#pragma once
/******
 * Purpose: This program is design to do mapping SOLiD or Illumina Reads,
 * The input is a set of short reads file and a Genome_Indes_TableQ
 * The mapping result is output to a file
 *  Author: Yangho Chen
 */
#include "LongReadsSet.h"
#include "ReadsMappingStats.h"
#include "MappingResult.h"
#include "Genome_Index_TableQ.h"
#include "ParameterList.h"
#include "AlignmentsQ.h"
#include "ReadInBitsSet.h"
#include "ReadInBits.h"
#include "ShortReadUtil.h"
#include "ColorSpaceRead.h"
#include "GenomeNTdata.h"
#include "chromosomeNTdata.h"
#include "stdafx.h"
#include <time.h>
#include <vector>
#ifdef WIN32
#include <direct.h>
#endif

// Use OpenMP is gcc version is later than 4.2
#ifdef __GNUC__
#ifdef __GNUC_PATCHLEVEL__
#define __GNUC_VERSION__ (__GNUC__ * 10000 \
                            + __GNUC_MINOR__ * 100 \
                            + __GNUC_PATCHLEVEL__)
#else
#define __GNUC_VERSION__ (__GNUC__ * 10000 \
                            + __GNUC_MINOR__ * 100)
# endif

#if __GNUC_VERSION__ >= 40200
#include <omp.h>
#endif
#else
#ifdef _MSC_VER
#if _MSC_VER > 1500
#include <omp.h>
#endif
#endif
#endif

// Macro for Parallelization with OpenMP
#ifndef _OPENMP
#define __OPENMP_FOR_PARALLEL__(SHARP_SIGN, openmp_flag) \
{LOG_INFO("Info %d: Use single CPU because OpenMP is not available.\n", CONFIG_LOG);}
#else
#define __OPENMP_FOR_PARALLEL__(SHARP_SIGN, openmp_flag) {\
    int numberOfCPUs = omp_get_num_procs();\
    std::cout << numberOfCPUs << " CPUs." << BLANK_LINE << "\n";\
    SHARP_SIGN##openmp_flag omp parallel for\
}
#endif

// For both counter for chromosome and tolerated substitution error
const int ALIGNMENT_RESULT_FILE_BUFFER_SIZE = 1000000;

/*
 * This class maps Illumina or SOliD short reads to reference genome.
 * It switches different ways to do read mapping, according to different options.
 */
class CReadsMapping: public CReadsMappingStats
{
public:
    CReadsMapping(void);
    CReadsMapping(const MappingOpts P);
    ~CReadsMapping(void);
    MappingOpts opt;
    char cOutputFormat;

protected:
    int printMapInfo(CReadInBitsSet& readsSet, int seedOpt);
    int setUpIO4Aligment(const char* Solexafile, const CGenome_Index_TableQ& table);
    int tearDownIO4Aligment(void);
    FileOutputBuffer*  AlignResult;
    FileOutputBuffer*  MissReads;
    time_t start, end;

    const static unsigned int CHECK_POINTS = 1000000;
    int dealMappedRead(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue);
    int dealMissedRead(bool bMapReadInColors, const char* tag, CReadInBits r, const char* qs = NULL);
    int dealMissedRead(CMappingResult& m);
private:
    void initialization(void);
    unsigned int checklongReadSetSize(CReadInBitsSet& firstHalfSet, CReadInBitsSet& SecondHalfSet);
};
// The following function fill different part of CMappingResul
inline
void getSingleMappingIndex(CGenomeNTdata& pgenomeNT, CAlignmentsQ& aQue, int mappingIndex, CMappingResult &m)
{
    m.uiDiff = aQue.asdiff[mappingIndex];
    m.MultipleMappedNo = aQue.load;
    m.strand = (mappingIndex >= (int)aQue.ForwardAlignmentLoad) ? '-' : '+' ;
    m.uiGlobalMappedPos = aQue.aiHitIndex[mappingIndex];
    m.uiRefId = pgenomeNT.genomeIndex2chrID(m.uiGlobalMappedPos);
    if (m.uiRefId >= pgenomeNT.iNo_of_chromosome) {
        ERR
    }
    m.uiPOS = pgenomeNT.genomeLocusID2chrIndex(m.uiGlobalMappedPos);
    unsigned int chrLength = pgenomeNT.paChromosomes[m.uiRefId]->iChromosome_size;
    ChrIndex2GeneName& geneVec = pgenomeNT.paChromosomes[m.uiRefId]->geneVec;
    if (geneVec.table.size() > 0) {
        CGene g = geneVec.query(m.uiPOS);
        m.uiRefLength = geneVec.getGeneLength(geneVec.getGeneId(m.uiPOS), chrLength);
        strcpy(m.RNAME, g.name.c_str());
        // g.startIndex is the tranlated index of m.uiPOS. Not the start index of gene
        m.uiPOS = g.startIndex;
    } else {
        sprintf(m.RNAME, "%d", m.uiRefId);
    }
}

void getReadQscores4Solexa(CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getSingleMappingSeqAndQ4SOLiD\
(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue, CMappingResult& m, bool samFormat);
void getSingleMappingSeq4Solexa(const CGenome_Index_TableQ& table, CMappingResult& m, bool samFormat);
void getLongMappingInfo(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue, bool samFormat,\
                        unsigned int mappingId, CMappingResult& m);
void getSingleMappingInfo(const CGenome_Index_TableQ& table, CAlignmentsQ& aQue,\
                          unsigned int mappingId, CMappingResult& m, bool samFormat);

inline void printSamHeader(FileOutputBuffer* AlignResult, vector<CGene>&refs, const char* RG, const char* CL)
{
    sprintf(AlignResult->caBufp, "@HD\tVN:0.1.5c\tSO:queryname\n");
    AlignResult->UpdateSize();
    for (vector<CGene>::iterator it = refs.begin(); it != refs.end(); it++ ) {
        sprintf(AlignResult->caBufp, "@SQ\tSN:%s\tLN:%u\n", it->name.c_str(), it->startIndex);
        AlignResult->UpdateSize(); // startIndex is actually the length of reference.
    }
    sprintf(AlignResult->caBufp, "%s\n@PG\tID:Clippers\tVN:0.0.8\tCL:\"%s\"\n", RG, CL);
    AlignResult->UpdateSize();
}

bool wrongIndex(const CReadInBitsSet& readsSet, const CGenome_Index_TableQ& table);

